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RELECTION  OF  X-RAYS  FROM  REPEATED  MULTILAYER  STRUCTURES 


1.  Introduction. 

Shortly  after  x-rays  were  discovered  nearly  100  years  ago,  it  was  observed  that  crystals  were  able  to  reflect 
x-rays,  but  only  in  certain  directions.  The  accepted  explanation  for  this  phenomenon  -  called  x-ray  diffraction  - 
postulates  that  die  crystal  consists  of  simple  but  infinitely  repeated  identical  ’cells*  of  atoms;  for  a  given  wavelength 
and  an  arbitrary  direction,  reflections  from  the  many  cells  interfere  with  each  other;  it  is  only  along  a  few  directions 
that  the  reflections  from  the  many  repeated  cells  cooperate  (reinforce)  to  provide  a  detectable  ’diffracted*  beam. 
The  intensity  along  each  diffracted  beam  depends  on  the  exact  position  of  each  atom  in  the  cell.  Two  uses, 
conceptually  each  other’s  opposite,  result  from  these  facts:  on  the  one  hand,  crystals  can  be  used  to  produce  desired 
reflections  of  x-rays;  on  the  other,  observed  reflections  can  be  used  to  determine  the  detailed  atomic  structure  of 
the  cell.  An  early,  but  authoritative,  description  of  these  effects  appears  in  ref.  [1]. 

It  will  be  plausible  that  the  phenomenological  description  and  theory  summarized  above  will  be  appropriate 
only  if  the  size  of  the  cells  and  the  wavelength  of  the  radiation  are  of  roughly  the  same  size:  on  the  one  hand,  very 
short  wavelengths  'will  not  notice’  the  correlation  between  atoms  in  a  cell,  while,  on  the  other,  very  long  ones 
won’t  even  notice  the  periodicity.  So,  if  you  are  interested  in  reflecting  long  wavelengths,  you  might  not  find 
natural  crystals  of  large  enough  cell  size;  you  might  have  to  construct  your  own.  This  is  the  rationale  of  the 
attempts,  in  the  last  few  decades,  of  constructing  repeating  ’multilayers*  as  artificial  crystals  [2].  These  are, 
generally  speaking,  a  thin  film  of  substance  A  of  precisely  known  thickness  (usually,  a  few  atomic  layers),  followed 
by  a  similar  layer  of  substance  B;  followed  by  more  identical  bilayers  ,  ABABABAB . 

The  diffraction  properties  of  repeating  multilayers,  usually  repeating  bilayers,  have  been  studied  at  the 
Naval  Research  Laboratory  (NRL)  both  experimentally  and  theoretically.  For  the  theoretical  work,  two  approaches 
have  been  used: 

1)  an  x-ray,  or  atomistic  approach,  and 

2)  an  electromagnetic,  or  homogeneous  approach. 

In  1),  a  long  and  thin  unit  cell  is  first  defined,  as  shown: 


P  . >x 

Both  a  and  b  atoms  are  in  "ordered*,  specified  positions.  This  long  and  thin  cell  is  repeated  an  infinite  number  of 
times  in  the  space  directions  x,y,  and  z  to  form  a  semi-infinite  slab.  The  reflective  properties  of  the  multilayer  is 
then  calculated  by  evaluating  the  structure  factor  and  other  procedures  well  known  to  workers  in  x-ray 
crystallography;  see  ref  [1].  Corrections  are  later  made  for  the  finiteness  of  the  multilayers,  the  vibrational  motion 
of  the  atoms,  and  the  absorptive  properties  of  the  layers.  The  fact  that  atoms  in  thin  layers  are  probably  in  random 
rather  than  ’ordered’  positions  is  taken  as  unimportant  for  long  enough  wavelengths.  For  details  of  this  atomistic 
approach,  refer  to  refs.  [3]  and  [4]. 

The  electromagnetic  approach,  2),  is  the  one  described  in  the  rest  of  this  report.  Each  layer  is  assumed 
to  be  homogeneous  (non-atomic),  the  incident  radiation  is  taken  to  be  an  electromagnetic  wave,  and  their  interaction 
described  by  Maxwell’s  equations.  It  is  thus  a  purely  "classical"  (i.e.,  a  non-quantum)  theory,  containing  neither 
atoms  nor  photons.  At  each  interface,  the  electromagnetic  wave  is  split  into  a  refracted  and  a  reflected  part  given 
by  Fresnel’s  laws  (  ref.  [4] )  (which,  of  course,  are  derivable  from  Maxwell’s  equations);  within  the  interior  of  any 
layer,  the  wave  is  attenuated  by  absorption.  We  can  thus  calculate,  successively,  the  properties  of  the 
electromagnetic  wave  after  any  number  of  bilayers  -  until  we  reach  the  thickness  of  the  specific  multilayer  we  want 
to  describe,  or  until  absorption  has  reduced  the  intensity  to  a  value  so  low  that  is  no  longer  interesting. 
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There  is  so  new  physics  in  this  method  of  calculating  the  effect  of  multilayers  on  electromagnetic  radiation 
(  see  refs.  [S],  [6],  [7] ) .  Rather,  what  is  done  here  is  the  development  of  a  formalism  and  a  computer  program 
conveniently  applicable  to  a  problem  of  continuing  interest  at  the  Naval  Research  Laboratory  [3]. 

Do  the  two  methods  agree?  If  not,  which  (if  either)  is  correct,  and  under  what  circumstances?  We  have 
only  partial  answers  to  these  very  reasonable  questions.  As  noted  above,  our  method  2)  is  more  likely  to  be  valid 
the  longer  the  wavelength  of  the  radiation;  for  short  waves  (’hard  x-rays’),  atomic  interactions  cannot  be  ignored 
or  averaged  over  and  a  version  of  method  1)  must  be  used.  On  the  other  hand,  method  2)  takes  more  reasonable 
account  of  several  physical  properties:  unlike  method  1),  it  does  not  have  to  assume  strict  periodicity  on  the  atomic 
level;  the  fact  of  absorption  enters  the  calculation  properly  ab  initio,  rather  than  as  a  correction  to  an  absorption- 
free  calculation;and  the  same  is  true  for  the  fact  that  the  number  of  bilayers  in  the  structure  under  consideration 
is  finite  rather  than  infinite.  In  addition,  method  2)  is  able  to  compute  the  reflectivity  at  any  wavelength  and  in 
particular  the  shape  of  any  reflection  line,  while  method  1)  gives  only  the  integrated  reflectivity. 

2,  Soft  X-rav  Reflection  via  Classical  Optics 

a)  a  single  layer 

As  sketched  in  figure  1,  consider  an  electromagnetic  wave  going  from  region  1  down  into  regions  2  and 
region  t  (mnemonic  for  future  use:  t  stands  for  "last").  Bom  and  Wolf  (  ref  [5];  referred  to  as  BW  )  show  that 
the  electromagnetic  field  at  point  h2  is  related  to  that  at  point  0  by  the  relation 


U(h2)  -  m2  .  1/(0) 


where  U(h)  is  defined  as  the  1x2  matrix 


“«■>-  (#>>] 


with  E  and  H  the  electric  and  magnetic  fields,  respectively, 
and  M2  is  the  2x2  matrix 

M  ~  \  C0S^2  -(i/p2)sin02l  (3a) 

2  "  [  ~(iP 2)^ 2  ] 


with 


P2  “  (n2/n2)  costfj 
/3j  —  (2x/jj/X,)n2  cosflj 
sin02  —  (;i,/n2)  0, 


(1) 


(2) 


Please  note  that  0  is  defined  as  the  angle  of  incidence  as  measured  with  respect  to  the  normal  to  the  material  surface, 
as  shown  in  figure  1  and  as  is  customary  in  literature  on  optics.  (  In  the  field  of  x-rays,  0  usually  denotes  the 
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complement  of  this  angle).  The  subscript  2  in  the  various  quantities  indicates  that  they  refer  to  material  *2*  .  The 
parameters  £,  n,  n  and  are  the  dielectric  constant,  the  magnetic  permeability,  and  the  index  of  refraction;  they  are 
related  by  iJir  „(. 

We  can  now  use  equ.  (1)  to  calculate  the  field  quantities  E  and  H  below  the  layer  in  terms  of  their  values 
above  it.  We  think  that  the  reader  will  find  it  reasonable  that  the  reflectivity  of  the  layer  of  material  *2*  which 
is  1>2  centimeters  thick  can  indeed  be  calculated  from  these  field  quantities.  The  details  •  a  little  lengthy  but 
straightforward  -  are  given  by  BW  in  their  equs.  (48)  and  (SI)  of  sec.  6.1.  The  polarization  of  the  incident  beam 
also  enters  into  the  calculation. 

b)  Bilayers  and  repeated  bilayers. 

The  virtue  of  the  matrix  formulation  (1)  is  that  the  effect  of  any  sequence  of  layers  can  now  be  written 
down  effortlessly:  suppose  that  the  layer  of  material  *2’  is  followed  by  a  layer,  h3  cm  thick,  of  material  *3",  as 
shown  in  figure  2  :  then  we  have 


U(h2+hJ  -  A/3  .  A/2  .  U(O) 


(4) 


where  M3  is  just  M2  with  subscripts  2  replaced  by  3;  and  for  a  sequence  of  N  bilayers  of  materials  2  and  3  we  have 
simply 


U[  N(h2+hJ  ]  -  (A/, .A//  .  U(0) 


(5a) 


as  illustrated  in  figure  3.  To  find  the  reflectivity  of  these  N  bilayers,  we  proceed  just  as  in  section  3a)  above, 
except  that,  in  solving  (5)  instead  of  (1),  we  must  use  the  more  complicated  (  but  still  2  by  2)  product  matrix 
(M3.Mj)n  instead  of  M2. 

For  convenience,  rewrite  (5a)  as 


l/( bottom )  —  A/£  .  U(top) 


(5b) 


where 


A/a  -  A/3  .  A/, 


(6) 


is  found  by  direct  multiplication,  and  is  written  down  by  BW  (their  p.67,  equ.  (86) )  and  also  in  our  appendix  1. 
It  can  be  written  in  the  form 


A/„ 


C 


c 

b) 


(7) 
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Next,  we  need,  for  equ.(5),  the  Nth  power  of  this  matrix.  We  do  this  by  diagonalizing  MB;  that  is,  we  find  the 
matrix  S  which  produces 


S'1  M  S- 


X,  0 
0  Xj 


(8) 


(where  we  have  for  simplicity  droped  the  subscripts  23), and  also  find  the  "eigenvalues’  X,  and  X2.  The  eigenvalues 
of  the  Nth  power  of  this  matrix,  needed  in  (5),  come  out  simply  to  be  X|N  and  X2N,  as  is  seen  from 


S'1  MNS  -  S'  MMM  ...  MMS 

-  S’1  M  (SS')  M  ( SS' )  ...  U  ( SS' )  MS 

-  (S-'MS)  ( S'MS)  ...  (S'MS) 

-  ( S'lMS)N 


which  with  the  use  of  equ.  (8)  becomes 


S'1  MN  5  - 


X,  0 

N 

xf  0 

0  Xj 

0  X? 

(9) 


or 


M" 


0 


X? 


S 


This  explains  why  the  number  N  of  multilayers  appears  in  appendix  1  in  such  a  simple  way  -  as  an  exponent  on  X, 
and  X2.  All  we  shall  need  are  the  explicit  values  of  X,.  X2,  and  S.  These  come  out  of  the  diagonalization  process 
and  are  given  in  appendix  2. 

We  now  proceed  just  as  in  the  last  paragraph  of  section  2a):  put  (9)  into  (5b)  to  get  the  field  quantities 
E  and  H  below  the  multiple  bilayer  from  their  values  above,  and  calculate  the  reflectivity  from  them. 

This  completes  the  essentials  of  the  calculation;  the  details  appear  in  the  appendices.  Appendix  1  covers 
the  mathematical  aspects  of  this  section.  Appendix  2  details  the  diagonalizing  transformation  of  a  2x2  matrix. 
Appendix  3  shows  how  the  existence  of  reflection  peaks  is  related  to  the  analytical  properties  of  the  eigenvalues  of 
the  characteristic  matrix.  Appendix  4  gives  two  explicit  expressions  for  the  reflectivity  r  (which  was  derived  in 
appendix  1),  and  notes  the  conditions  under  which  either  is  preferable  for  computations.  Appendix  5  relates  the 
data  for  each  atomic  constituent  of  a  layer  to  the  gross  properties  (i.e.  the  index  of  refraction)  of  the  layer. 
Appendix  6  reconciles  a  ..jiational  difference  between  two  references.  Appendix  7  is  a  printout  of  the  computer 
program  that  is  described  verbally  in  section  5. 

3.  Integrated  reflectivity. 

The  preceding  section,  together  with  the  details  in  the  appendices,  has  allowed  us  to  compute  the  reflectivity 
of  a  repeated  bilayer  at  a  specified  incident  wavelength.  The  obvious  next  step  would  seem  to  be  the  repetition 
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of  the  calculation  for  a  different  but  nearby  angle,  until  the  entire  range  of  tbeta  from  0  to  pi/2  has  been  covered. 
This  is  indeed  done  in  our  computer  program. 

A  well-known  experimental  fact  in  x-ray  theory  and  practice  is  the  existence  of  *  reflection  lines*;  that  is, 
strong  reflection  over  one  (or  more)  very  small  regions  of  wavelength,  with  (near-)zero  reflection  between  them 
(  see  ref.[l]  ).  For  a  sufficiently  large  number  N  of  bilayers,  our  calculations  verify  this  state  of  affairs.  Naturally 
one  then  wants  to  know  the  ’integrated  reflectivity*  of  one  line,  and  this  quantity  can  be  obtained  by  summing  (or 
'numerically  integrating’)  over  the  wavelengths  contained  in  one  line. 

We  have  of  course  incorporated  this  into  our  computer  program,  but  want  to  warn  that  inherent 
arbitrariness  remains  in  the  choice  of  the  limits  in  the  integration.  The  method  we  have  adopted  is  to  put  the  lower 
limit  of  integration  1/3  of  the  way  to  the  line  to  the  left,  and  the  upper  limit  1/3  of  the  way  to  the  line  to  the  right. 
This  will  be  fine  if  each  line  is  indeed  'sharp*  and  if  the  reflection  is  very  close  to  0  in  between;  but  exceptions 
to  this  rule  will  be  buried  beyond  recognition  by  this  choice  for  limits  of  integration.  Caution  is  advised. 

4.  Surface  Roughness. 

Experimental  observation  of  lines  that  are  broader  or  weaker  than  predicted  by  theory  have  been  plausibly 
attributed  to  'surface  roughness*;  see  ref.  [5].  Our  theoretical  model  describes  layers  with  two  properties: 

1)  boundaries  are  perfect  planes,  and 

2)  bilayers  are  repeated  with  perfect  periodicity. 

Neither  of  these  conditions  is  likely  to  be  fully  attained  in  the  real  world;  it  is  plausible  to  attribute  the  deviations 
of  experimental  data  from  theory  to  surface  roughness. 

Can  we  put  surface  roughness  into  our  model  while  maintaining  mathematically  essential  properties  1) 
and  2)?  What  we  have  done  is  to  replace  the  two  layers  consisting  of  pure  material  A  and  pure  material  B  by  eight 
layers  of  the  same  total  thickness;  layer  2  is  pure  A,  layer  6  is  pure  B,  and  layers  3,4,5  have  intermediate  index 
of  refraction,  as  do  layers  7,8,9.  The  procedure  is  justified  in  greater  detail  in  ref.  [8]. 

5.  Description  of  Computer  program. 

A  version  of  the  FORTRAN  computer  program  in  use,  called  LA8WC.FOR,  is  presented  in  Appendix  7. 
The  line  numbers  on  the  far  left  appear  for  the  reader’s  convenience  only  and  are  ignored  by  the  computer;  the  set 
of  numbers  appearing  to  the  right  of  the  line  numbers  are  FORTRAN  statement  numbers.  The  program  displayed 
in  Appendix  7  describes  a  repeated  structure  of  8  layers,  numbered  from  2  to  9;  layer  2  consists  of  W  (tungsten) 
and  layer  6  of  C  (carbon);  the  other  layers  have  an  intermediate  composition,  consisting  of  tungsten  and  carbon  ions 
in  the  ratios  3:1,  1:1,  and  1:3.  The  program  can  be  applied  to  other  multilayer  materials  by  changing  a  small 
number  of  lines  between  lines  50  and  166. 

For  easier  readability,  we  always  attach  a  statement  number  n  which  ends  in  0  to  the  first  statement  in  each 
"do  loop*,  and  statement  number  n  + 1  to  the  last  statement  in  that  loop,  e.g.  thus: 

800  do  801  kk—  1,  9 


801  continue. 

The  main  part  of  the  calculation  are  the  nested  lops  which  start  at  line  number  245: 
700  sums  over  the  "orders"  of  the  reflection  peaks 

200  sums  over  the  two  polarizations  (TE,  TM) 

300  sums  over  the  angles  of  incidence  around  a  peak 


5 


400  sums  over  the  layers  in  one  cell 

401 
301 

201 

701 

The  main  printout,  which  is  the  integrated  reflectivity  of  one  peak,  is  elicited  by  statement  #  50  (  at  line  #  388). 
However,  the  (non-integrated)  reflectivity  at  each  angle  of  incidence  can  also  be  printed  out,  by  removing  the  c  (= 
"comment*)  from  statement  66  (at  line  number  382). 

The  244  lines  that  precede  this  main  calculation  will  not  be  described  in  much  detail  here,  since  they  are 
comparatively  simple  structurally,  proceeding  in  a  linear  fashion  without  much  nesting  or  interrelationships.  Lines 
1  through  33  are  explanatory  comments;  34  through  49  are  FORTRAN  declarations  of  variables.  Basically  what 
is  done  in  lines  SO  through  166  is  the  insertion  of  the  data  specific  to  the  chemical  species  involved;  they  culminate 
in  the  calculation  of  the  index  of  refraction  (called  ninx)  in  lines  173-  195,  which  is  the  quantity  used  in  the  four 
main  loops. 

A  possible  point  of  confusion-  viz.  a  change  in  numbering  of  the  layers  -  is  described  in  lines  205,  206. 
The  reason  for  this  is  historical,  not  logical:  this  program  for  the  8  layers  per  cell  was  constructed  from  an  earlier 
one  for  4  layers  per  cell.  The  program  could,  of  course,  be  rewritten  to  give  each  layer  the  correct  number  in  the 
first  place. 

We  have  in  many  cases  provided  two  ways  of  inserting  data,  at  the  option  of  the  operator:  from  the 
keyboard,  or  by  modifying  a  statement  in  the  program.  For  example,  in  lines  95-100,  densities  are  inserted  by 
statements  in  the  program.  To  change  to  insertion  from  the  keyboard,  remove  the  ”!’  from  line  97,  and  insert  a 
"!*  in  lines  98  and  99. 

To  run  the  program,  the  command,  to  be  entered  from  the  "$"  prompt,  is 
©exnoop  la8wc 

This  calls  a  short  command  file,  EXNOOP.COM  ,  which  provides  the  usual  FORTRAN,  link,  and  run  commands 
without  use  of  the  optimizer.  The  reason  for  the  exclusion  of  the  optimizer  is  that  it  often  gives  wrong  answers. 
This  is  a  problem  of  the  computer  we  are  using  that  will,  we  are  told,  be  fixed  in  the  future. 

Appendix  1.  Mathematical  details. 

This  appendix  covers  the  same  ground  as  sec.  2,  but  does  so  in  detailed  mathematical  language,  and  in 
notation  similar  to  that  used  in  the  computer  program  (Appendix  7),  with  minimal  verbal  explanation. 

The  results  depend  on  the  polarization  of  the  incoming  beam.  The  formulas  below  refer  to  TE  polarization. 
TE  (=  "transverse  electric"]  means  that  the  Electric  field  vector  E  is  perpendicular  ("Transverse")  to  the  plane  of 
incidence,  i.e.  to  the  plane  of  the  paper  in  figure  1.  TM  is  analogously  defined  for  the  Magnetic  field  H.  At  the 
end  of  this  appendix,  we  explain  how  the  TM  formulas  can  be  easily  obtained  from  the  TE  ones.  To  describe  an 
unpolarized  beam,  the  TE  and  TM  results  should  be  averaged  in  the  end. 

The  index  i  refers  to  the  materials  involved.  Orginally,  i=2  and  3  referred  to  the  bilayer,  which  is 
repeated  N  times;  i-  1  and  (  are  not  repeated;  they  refer  to  the  materials  bounding  the  repeated  bilayers.  At  the 
end  of  this  appendix,  we  explain  how  the  formulas  must  be  modified  when  the  repeated  structure  consists  of  more 
than  2  layers. 

Input  data  are:  angle  of  incidence  6,,  incident  wavelength  X0,complex  indices  of  refraction  n,,  layer 
thickness  h2,  h3,  number  of  bilayers  N  (called  layno).  Note  that  generally  nj=  (ej/*j)'/i.  where  t  and  n  are  the 
dielectric  constant  and  the  magnetic  permeability.  We  first  compute  h=  h2  +  h3  and  and  then  the  sin  6,  =  (n./nj 
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sin  6 ,.  This  is  Snell's  law;  note,  however,  that,  as  nj  is  generally  complex,  sin  0,  is  also,  contravening  the  simple 
interpretation  of  6  as  an  "angle*.  Further: 


Pi  *  (nj/pi)  cosfli 

ft  =  (2t/>0  hjiij  cosflj 

a  =  cos02  cos/?3  -  (pj/pJ  sm02  sin03 

b  =  cos02  cos0y  -  (pj/pj)  sin02  sin03 

c  =  -i[(l/p3)  cosj32  sin/J3  +  (I/P2)  sinj32  cos/?3] 

d  =  -i[  p3  cosffi  sin/?3  +  p2  sin/?2  cos/J3] 

a  =  b-a  =  [(p3/p2)  -  (P2/P3)]  sin/32  sinb/33 

T  =  [a2  +  4cdf 

f  =  a  +  T 

X,  =  (lA)  [a  +  b  +  T] 

X2  =  ('A)  [a  +  b  -  T] 

D  =  2  fT 

M„  =  (4cdX2N  +  f  X2N)/D 

M22  =  (4cdX2N  +  f2  X,N)/D 

M,2  =  2cf(X,N -X2n)/D 

Mj,  =  2df  (X,N  -  X2n)/D 

P  =  (M,|  +  M,2  p()  p, 

Q  =  (M2i  +  mi3)  pf 

r  =  (P-Q)/(P  +  Q)  =  reflection  coefficient 

R  =  |  r  | 2  =  reflectivity 


The  reader  of  sec.  2  and  appendix  2  will  recognize  the 


matrix 


as  giving  the  values  of  E  and  H  below  a  bilayer  in  terms  of  the  values  above  it;  the  matrix  arises 

from  mdltiplyfag  the  characteristic  matrices  of  layer  2  and  layer  3.  This  tells  us  how  to  modify  the  calculations 
when  more  than  two  —  e.g.  three  —  different  layers  are  present:  simply  replace  equ.  (6)  by 
M234=  M4  M3  Mj,  and  so  on  for  any  number  of  different  layers. 


The  equations  below  a,  b,  c,  d  extract  the  reflectivity  from  the  field  quantities  E  and  H  above  and  below 
the  N  multilayers,  as  qualitatively  explained  in  sec.  2. 

We  noted  that  the  above  formulas  apply  to  TE  polarized  radiation.  To  obtain  equivalent  results  for  TM, 
only  one  change  must  be  made:  for  all  i,  replace  ^  —  (njn) cos(0()  by  pt  ( ji(/n()cos(fl ,).  See  ref.  [9]. 


Appendix  2.  Diaeonalization  of  a  2  bv  2  matrix. 

The  process  of  diagonalizing  an  n  by  n  matrix  is  well  known,  in  the  sense  that  is  described  in  many  text 
books  (  e.g.  refs  [10],  [11])  and  carried  out  in  several  published  computer  programs.  Numerical  methods  must 
generally  be  used,  either  from  the  beginning  or,  at  any  rate,  before  the  end.  However,  the  process  becomes  much 
simpler  for  2  x  2  matrices:  the  secular  equation  is  then  quadratic,  and  the  eigenvalues  are  thus  explicitly  obtainable 
in  terms  of  radicals,  as  are  all  other  quantities.  Hence,  everything  can  be  done  analytically.  Not  having  found  the 
simple  2x2  case  described  explicitly  in  the  literature,  we  summarize  the  results  here.  We  do  not  derive  them, 
since  the  reader  can  easily  verify  that  the  S  given  below  does  indeed  diagonalize  M  according  to  equ.  (8). 

The  matrix  of  interest  is  given  by  equ.  (6), 

# 


Its  inverse  is 


M'x 


/( ab-cd ) 


The  eigenvalues  defined  by  equ.  (8)  are 


X,  “  ( a+b+T)l2  and 
X,  “  (a+b-T)/2,  where 

T  “  UI*4a#]wl  and 
a  —  b-a. 


The  transformation  (8)  is  accomplished  by  the  matrix 


-/ 

2d 


and  its  inverse  is 


where 


We  begin  by  writing  down  the  matrix  describing  one  single  layer  of  index  n  according  to  equ.  (3)  of 


section  2: 


w  _  /  cosy  -(Up)  siny  \ 

\  -(ip)siny  cosy  / 


(A3-1) 


with 

y  —  (2* h/\J  cos6c 
p  —  (nip)  cos 6t 

We  note  that  this  is  a  "unimodular"  matrix  (which  means  that  its  determinant  =  1).  The  matrix  describing  any 
number  of  layers,  possibly  different  ones,  is  therefore  also  unimodular,  being  a  product  of  matrices  of  the  above 
form;  and  this  is  also  true  for  any  similarity  transform  H  M  H'1  of  any  such  product  matrix,  since  the  determinant 
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of  H  and  H'1  are  each  other’s  reciprocal;  in  particular  it  is  true  for  the  diagonal  matrix  whose  elements  are  the 
eigenvalues  X,  and  X2.  That  is,  the  product  of  the  two  eigenvalues  obeys 


\  \  -  1. 


(A3-2) 


This  relation  has  a  different  meaning  depending  on  whether  the  eigenvalues  are  real,  pure  imaginary,  or  complex. 

Real  eigenvalues.  If  one  X  is  real,  then  (A3-2)  implies  that  the  other  is  also,  and  is  the  first  one’s 
reciprocal;  and  one  X  is  greater  than  1,  the  other  smaller.  E.g.,  if  X  =  8,  then  X?  =  0.12S. 

Pure  imaginary  eigenvalues.  Here  (A3-2)  implies  that  if  one  is  pure  imaginary,  the  second  one  is  also. 
If  we  write  X,  =  iu,,  X2  =  iuj  (where  the  u  s  are  real),  then  (A3-2)  implies  that  X2  =  -i/u,.  E.g.  if  X,  =  10  i,  then 
Xj  -  -i/10. 

Complex  eigenvalues.  These  can  be  written  X,  =  r,*exp(ia),  X2  =  r2*exp(ib)  (where  a,  b,  and  the  r’s  are 
real).  To  satisfy  (A3-2),  we  must  have  b=  -a,  so  we  have 

X,  —  r,*exp(ia),  \  —  r2  *  exp(-/a).  (A3-3) 

Now  recall  that  the  trace  of  a  matrix,  like  the  determinant,  is  invariant  under  a  similarity  transformation.  We  can 
therefore  equate  the  trace  after  the  transformation  to  the  trace  before, 

X,  +  Xj  -  TriM) 


or,  using  (A3-3)  for  the  left  and  (A3-1)  for  the  right. 


r,  exp (ia)  +  rz  exp (-id)  ■■  2  cosy. 


(A3-4) 


If  y  is  real,  then  cos  y  is  also,  and  it  follows  that  rl  and  r2  are  equal,  since  otherwise  laml  +  lam2  would  not  be 
real.  So  in  that  case  we  have 


X,  —  r  exp  (ia)  and  X,  —  r  exp(-ia),  (A3-5) 

i.e.  they  are  each  others  conjugates. 

What  is  the  physical  meaning  of  a  real  trace?  From  (A3-1)  we  see  that  Tr(M)  is  real  iff  n=  (€  is  real, 
i.e.  if  there  is  no  absorption  in  that  layer.  The  matrix  M  in  that  case  is  of  the  form 

real  imag 
imag  real 

Now  multiply  that  matrix  by  another  of  the  same  structure  (Physically:  follow  tbat  layer  with  another  non-absorbing 
one).  Then  the  product  matrix  has  the  same  form  also.  (This  is  perhaps  not  obvious,  but  easily  seen  by  carrying 
out  the  multiplication).  That  form  will,  by  the  same  reasoning,  persist  through  any  number  of  multiplications 
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(physically:  through  any  number  of  non-absorbing  layers).  We  conclude  that  the  simple  relation  (A3-5)  holds  for 
a  multilayer  consisting  wholly  of  non-absorbing  layers,  while  in  presence  of  absorption  only  the  more  general 
relation  (A3-3)  applies. 

Now  what  is  the  purpose  of  all  this?  We  want  to  find  the  behavior  of  the  multilayer  calculation  as  a 
function  of  N,  the  number  of  double  layers.  N  appears  in  the  present  calculation  only  in  the  form  X,N  and  X2N. 
The  limit  of  large  N  is  particularly  important,  as  is  stability  :  for  physical  reasons,  a  constant  value  should  be 
approached  when  N  gets  large. 

Neither  complex  eigenvalues,  nor  pure  imaginary  ones,  provide  that  convergence:  X,N  has  the  form 
exp(iaN),  which  changes  substantially  with  every  unit  increase  of  N;  XjN  behaves  in  the  same  way.  The  same  is 
true  for  pure  imaginary  eigenvalues.  On  the  other  hand,  for  real  eigenvalues,  the  greater  of  the  two  will,  when 
raised  to  the  power  N,  be  much  larger  than  the  smaller  raised  to  the  same  power.  The  expression  for  the 
reflectivity  derived  in  App.  1  contains  the  eigenvalues  in  both  the  numerator  and  the  denominator,  and  the 
reflectivity  r  then  becomes  independent  of  both  eigenvalues,  hence  also  independent  of  N;  this  is  perhaps  best  seen 
from  eqs.  (A4-3).  In  absence  of  absorption,  that  limiting  value  is  1,  the  largest  possible  value.  (This  requires  a 
short  calculation.)  In  the  special  case  that  neither  eigenvalue  is  larger  than  the  other,  viz.  that  they  are  equal,  r 
vanishes. 

We  conclude  that  the  large  r’s  will  arise  for  values  of  the  angle  of  incidence  0,  for  which  the  X’s  are  real. 
To  find  these  regions  of  0„  we  wrote  a  preliminary  program  called  EIGV.FOR,  which  prints  the  eigenvalues  as  a 
function  of  thetal.  This  we  followed  by  the  final  program  LAYER.FOR.  Together,  they  verify  the  above 
reasoning:  large  values  of  r  appear  only  where  the  eigenvalues  are  real. 

We  would  also  expect  that  increasing  N  would  increase  the  reflectivity,  and  would  sharpen  the  width  of 
the  line  (i.e.  decrease  the  range  of  thetal  for  which  reflection  is  substantial).  In  absence  of  absorption,  the 
reflectivity  approaches  1  as  N  approaches  infinity;  if  absorption  is  present,  a  finite  value  for  the  absorption  is 
approached  for  N  large  enough  to  absorb  essentially  all  the  incoming  radiation,  and  no  further  change  in  the 
reflectivity  should  result  from  increasing  N  further.  We  have  verified  these  features  by  appropriate  model 
calculations. 

Appendix  4.  Two  expressions  for  the  reflectivity  suitable  for  computation. 

As  the  reader  can  check,  the  expression  for  r  arrived  at  in  appendix  1  can  be  written  as 


(A4-1) 


where 


-2d  *  pj 

*2  “f  -  2P*C  (A4-2) 

K  ~2cPl±f 

K  -fPi±  2d 


However,  in  most  situations  of  interest,  one  of  the  two  terms  in  the  numerator  of  (A4-1)  will  be  much  larger  than 
the  other;  and  the  same  for  the  denominator.  This  follows  from  two  facts  derived  in  appendix  3:  high  values  of 
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r  arise  only  when  the  X’s  are  real;  and  the  two  X’s  are  each  other’s  reciprocal.  It  follows  that  one  X  will  be  larger 
than  1  and  the  other  smaller;  for  the  moment,  let  us  call  them  X,  and  X,.  If  N  is  a  large  number  (as  often  it  is), 
X,N  will  be  much  larger  then  X.N,  proving  the  verbal  statement  following  equ.  (A4-2).  It  is  therefore  a  good  idea 
to  divide  both  the  numerator  and  the  denominator  by  X,N.  We  obtain 


K  +  xw  (*2W4. 


(A4-3a) 


suitable  when  |  X,  |  >  |  X:  |  and 


*4.  *  WWl- 


(4-3b) 


suitable  when  |  X2 1  >  |  X,  | .  Note  that  both  (A4-3a)  and  (A4-3b)  are  exact  in  all  cases;  but  in  most  cases  only 
one  of  them  will  allow  the  computer  to  proceed  without  complaining  about  "overflow"  and  "underflow". 

Appendix  5.  Index  of  refraction  n  and  atomic  scattering  factor  f. 

Our  calculation  requires  the  index  of  refraction  for  each  of  the  atomic  species  involved  and  at  the  incident 
wavelength  of  interest.  We  use  the  relationship  given  by  James,  ref  [12]: 


n  -  1-6  -  1  -  (JVX2  eVlr  mc^O) 


(A5-1) 


Here  N  is  the  number  of  atoms  per  unit  volume,  X  the  incident  wavelength,  e  and  m  the  charge  and  mass  of  the 
electron,  and  c  the  speed  of  light;  f  is  the  atomic  scattering  factor,  and  the  argument  (0)  denotes  grazing  incidence. 

f  can  be  obtained  in  two  ways.  The  preferred  way  uses  the  tables  of  Henke  et  al.,  ref.  [13],  where  two 
quantities,  f,  and  f2,  are  tabulated  for  94  atomic  species  and  many  wavelengths.  They  are  the  real  and  imaginary 
parts  of  f(0).  That  is, 


where 

„d 

Z  is  the  atomic  number,  and  f  +  if”  is  commonly  called  the  "anomalous  scattering  factor". 

For  wavelengths  not  covered  by  Henke,  the  scattering  factors  are  taken  from  Cromer  [14],  who  gives  f 
and  f  in  his  equs  (8)  and  (9).  However,  a  numerical  integration  is  required  to  evaluate  equ  (7)  (in  contrast  to 
Henke’s  data,  which  only  need  to  be  looked  up). 

As  a  practical  matter,  a  user-friendly  computer  program  due  to  D.B.  Brown,  called  XTALR.COM  exists, 
ref.  [3],  which  in  its  preliminary  stages  finds  f  in  just  the  manner  described  above  -  i.e.  by  using  Henke’s  numbers 
when  available,  and  computing  them  according  to  Cromer  when  not.  Considerable  time  may  be  saved  by  utilizing 
that  program. 
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Finally,  a  word  about  the  index  of  refraction  of  a  layer  containing  several  different  atomic  species.  Since 
each  atom  contributes  additively,  relation  (A5-1)  can  be  used  with  two  modifications:  substitute  E  fj(0)  for  f(0), 
where  the  index  i  sums  over  all  the  atoms  in  a  "cell”;  and  reinterpret  N  as  the  number  of  "cells’  per  unit  volume. 

Appendix  6.  Sign  change  for  f. 

Our  analysis  (  Sec.2  )  is  based  on  BW,  ref.  [4],  but  we  also  use  formulas  from  J  ,  ref.  [12],  and  data  from 
H,  ref.  [13].  We  must  therefore  make  sure  that  the  notation  in  these  three  papers  is  consistent;  or  to  make 
appropriate  changes  where  it  is  not. 

We  have  found  one  inconsistency:  BW  and  H  write 
Ref  exp{+inkx)  ] 
for  a  plane  wave,  while  J  writes 

Ref  exp{-inkx)  ]. 

Since  J  is  the  one  who  is  out  of  step,  the  easiest  way  to  make  our  calculation  consistent  is  to  modify  J’s  formulas 
whenever  they  are  used:  we  must  replace  his  Im(n)  by  -Im(n).  A  detailed  validation  of  this  procedure  is  given  in 
ref.  [15]. 


Appendix  7.  Printout  of  computer  program. 

We  present  here  a  printout  of  a  computer  program  called  LA8WC.FOR.  The  many  comments  should  help 
its  readability;  however,  the  summary  in  sec.  5  of  this  report  should  prove  more  coherent. 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 


c  Program  LA8WC.FOR 

c  H.  B.  Rosenstock,  1990-1991 


c  program  with  eight  layers-  W,  C,  and  6  mixed  ones 
c  cgs  units  unless  otherwise  stated, 

c  General  notation  according  to  Bom  and  Wolf, 
c  mu,  eps,  ninx  =mag.  permeability ,diel.  constant, index  of  refr. 
c  theta(i)  are  angles  of  beam  with  the  normal  (not  with  the  plane!)  in  material  i. 
c  thel(lOOO)  are  1000  values  of  theta(l)  above 
c  h(2),  h(3)  are  thickness  of  materials  2  and  3, 

c  wavelengthO,  k0=  wavelength,  wave-number 

c  ha,  wavelengthangO,  kangO  =  same  in  Angstroms 
c  subno=  #  of  components(single  layers)  in  one  multilayer 

c  (not  counting  vacuum  at  both  ends) 

c  layno=  N=  0  of  multilayers;  Ndens=  #  of  atoms/vol;  Navo=  Avogadro’s  no.; 
c  Zat,Aat=  at.  charge, mass;  Zat  is  not  used  in  the  calculation 

c  proper,  but  may  be  needed  to  call  the  Henke  tables, 

c  ii=  sqrt(-l) 

c  300-301  theta-loop 

c  400-401  materials  loop  (inside  theta-loop) 

c  350-310  -  see  BW  sec.  6.1,  HBR  notes,  also,p(4),  beta(4) 
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27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 


c  500-loop  provides  printout  for  300  loop,  printing  maxima  and 
c  minima  of  reflectivity  only. 


implicit  real(a-z) 

complex  ii,  c,  d,  discri,  f,  laml,  lam2 
complex  refl 

complex  ninx,  sinthe,  costhe,  beta,  p,  theta,  decre 
complex  cb2,  cb3,  sb2,  sb3,  p2,p3,  a,  b,  delta 
complex  psil,  psi2,  psi3pl,  psi3mi,  psi4pl,  psi4mi 
complex  termpl,  termmi,  top,  bot,  fstar,  fsum,  tran,  prod 
panuneterfmax =26) 

dimension  mu(max),  ninx(max),  Asum(max),  fsum(max) 
dimension  h(max),  ha(max) 

dimension  theta(max),  sinthefmax),  costhe(max),  rho(max) 

dimension  Aat(max,10),  Zat(max,10) 

dimension  multy(max,10),  Ndens(max),  p(max),  beta(max) 

dimension  fl(max,10),  f2(max,10),  fstar(max,10),  thlpeak(0:10) 

dimension  thldiff(0:10) 

dimension  tran(2:max,  2,2),  prod(2:max,  2,2) 

ii=  (0,1) 
pi=  3.1415926 
Navo=  6.02e23  JAvogadro 
eel=  4.803e-10  !el.  charge 
mel=  9.l09e-28  !el.  mass 
clight=  3.00el0  ! speed 

0  type  1 

I  format  (’  write  a  label’) 
accept* 

c  layno=  1 

3  print  4 

4  format  (’  type  the  number  of  distinct  layers;  and 

2  of  multilayers  ’) 

!  accept*,  subno,  layno 
subno=  8 
layno =  600 

10  print  11,  subno,  layno 

I I  format  (’  number  layers  in  one  multilayer; 

2  and  of  multilayers  N  =’,  2f6.0  ) 

tyP«  2 

dim=  26  !not  to  exceed  stated  dimensionality  of  ninx, mu,  etc 
120  do  121  i  =  1,  dim  !  initial  settings;  change  them  later 

ninx(i)=  1  !  i  sums  over  the  9  layers  (1,9=  vacuum) 
mu(i)  =  1 
ha(i)  =  0 

130  do  131  kk=  1,  9  !  kk  sums  over  atoms  in  layer  i 

fl(i,kk)=  0 
12(i,kk)=  0 
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79 

Aat(i,kk)=  0 

80 

Zat(i,kk)=  0 

81 

multy(i,kk)=  0 

82 

131 

continue 

83 

121 

continue 

84 

5 

type  6 

85 

6 

fonnat(  ’  type  atomic  numbers  Zat(2,l),  (4,1),  (4, 2);  then 

86 

2  same  for  atomic  masses  Aat  ’) 

87 

» 

accept*,  Zat(2),  Zat(4),  Aat(2),  Aat(4) 

88 

Zat(2,l)=  74 

89 

Zat(4,l)=  6 

90 

! 

Zat(4,2)=  7 

91 

Aat(2,l)=  183.9 

92 

Aat(4,l)=  12.01 

93 

i 

Aat(4,2)=  14.01 

94 

type*  Zat(2,l),Zat(4,l),Zat(4,2),Aat(2,l),Aat(4,l),Aat(4, 

95 

7 

type  8 

96 

fonnat(  ’  type  densities  rho(2)  and  rho(4)  in  grams/cc’) 

97 

i 

accept*,  rho(2),  rho{4) 

98 

rho(2)=  19.3 

99 

rho(4)=  2.00 

100 

type*,  rho(2),  rho(4) 

101 

print  2 

102 

103 

type  2 

104 

37 

type  38 

105 

38 

format  (’  for  layer  2,  type  fl,  and  f2  ’) 

106 

!  accept*.  fl(2,l),  f2(2,l) 

107 

fl(2,l)=  43.25 

108 

f2(2,l)=  11.54 

109 

type*,  fl(2,l),  f2(2,l) 

110 

47 

type  48 

111 

48 

format  (’  for  layer  3,  type  the  fls,  and  the  f2s  ’) 

112 

!  accept*,  fl(4,l)  f2(4,l),  fl(4,2),  f2(4,2) 

113 

fl(4,l)=  6.24 

114 

f2(4,l)=  .305 

115 

fl(4,2)=  0.0  <7.12 

116 

f2(4,2)=  0.0  11.96 

117 

type*,  fl(4,l),  f2(4,l),  fl(4,2),  f2(4,2) 

118 

119 

type  2 

120 

multy(2,l)=  1 

121 

multy(4, 1)=  1 

122 

multy(4,2)=  0 

123 

27 

type  28,  multy(2,l),  multy(4,l),  multy(4,2) 

124 

28 

format(’  multiplicities’,  3f5.0 ) 

125 

type  2 

126 

2 

format  (’  ’) 

127 

128 

21 

type  22 

129 

22 

formate  type  incident  wavelength  in  Angstroms’) 

130 

!  accept*,  wavelengthangO 
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131 

132 

133 

134 

135 

30 

136 

31 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

160 

157 

161 

158 

159 

160 

110 

161 

i 

162 

163 

164! 

165 

Ill 

166 

167 

150 

168 

151 

169 

170 

171 

40 

172 

41 

173 

100 

174 

175 

176 

177 

178 

179 

180 

181 

800 

182 

wavelengthangO=  8.34 

lcangO=  2  *pi /wavelengthangO 

wavelengthO=  wavelengthangO*  le-8 

kO=  kangO*  le8 

type  31  ,  wavelengthangO,  kangO 

format  (’  wavelengthangO,  kangO  =  ’  2f6.2,  *  Angstroms,  M’) 
print  2 

!  ha(3)=  0.00 
!  ha(5)=  ha(3) 

!  ha(2)=  7.672  -  ha(3) 

!  ha(4)=  19.728  -  ha(3) 

!  ha(6)=  ha(2) 

!  ha(7)=  ha(3) 

!  ha(8)=  ha(4) 

!  ha(9)=  ha(5) 

ha345=  3.  !  total  thickness  of  layers  3  +  4  +  5 

ha(2)=  7.672  -  ha345 

ha(6)=  19.728  -  ha345 

ha(3)=  ha345/  3. 

ha(4)=  ha345/  3. 

ha(5)=  ha345/  3. 

ha(7)=  ha345/  3. 

ha(8)=  ha345/  3. 

ha(9)=  ha345/  3. 

type  161,  ha345 

format  (’  total  thickness  of  three  mixed  layers  is’,f5.2  ) 


hsum=  0 
do  111  i=  l.dim 

type*,  i,  ha(i) 
h(i)=  ha(i)*  le-8 
hsum=  hsum+  h(i) 
type*,  hsum 

continue 
hh=  hsum 
type  151 

format  (’  pause,  then  enter’) 
accept* 

print  41 

format  (’  i  mu  n  ha’) 

do  101  i=  2,  subno+1 

!  compute  index  of  refr.  acc.  to  R.W.James,  equ.(2.61) 
etanum=  wavelength0**2*  eel**2*  Navo 
etaden=  2*  pi*  mel*  clight**2 
eta=  etanum/  etaden 
Asum(i)=  0 
fsum(i)=  0 
decre=  0 
do  801  kk=l,9 

fstar(i,kk)=  fl(i,kk)-  ii*  f2(i,kk) 
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183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 


!  see  HBR  llJul89.rep  for  source  of  minus-sign  above 
fsum(i)=  fsum(i)+  multy(i.kk)*  fstar(i,kk) 

Asum(i)=  Asum(i)+  multy(i,kk)*  Aat(i,kk) 

801  continue 

!  next  three  lines  for  “vacant*  layer  (all  Aat  zero) 

805  if  (Asum(i)  .gt.  .1)  goto  806 

!  set  ninx(3)  and  (5)  equal  to  1 
dec  re  =  0 
goto  807 

806  decre=  eta*  rho(i)*  fsum(i)/  Asum(i) 

807  ninx(i)=  1-  decre 

101  continue 
type  2 

!  change  ninx(3)  and  (5)  equal  to  average  of  the  thick  ones 
!  ninx(3)=  (  ninx(2)+ninx(4) )/  2 

!  ninx(5)=  (  ninx(2)+ninx(4) )/  2 
!  ninx(6)=  ninx(2) 

!  ninx(7)=  ninx(3) 

!  ninx(8)=  ninx(4) 

!  ninx(9)=  ninx(5) 

!  above,  ninx(2)  describes  W,  and  ninx(4)  describes  C 
!  below,  ninx(2)  describes  W,  and  ninx(6)  describes  C 

ninx(6)=  ninx(4)  !  the  C  layer 

ninx(3)=  .75*ninx(2)  +  .25*  ninx(6) 
ninx(4)=  .50*  ninx(2)  +  .50*  ninx(6) 
ninx(5)=  .25*  ninx(2)+  .75*ninx(6) 
ninx(7)=  ninx(5) 
ninx(8)=  ninx(4) 
ninx(9)=  ninx(3) 

1100  do  1101  i=2,  subno+1 

!  print  the  new  values  for  all  i 
type  103,  i,  mu(i),  ninx(i),  ha(i) 

103  fonnat(  f4.0,  fl0.6,  ’  ’,  2fl0.6,  fl2.6) 

1101  continue 
print  2 

c  locate  reflection  peaks  from  Bragg’s  law;  and  the  distance 
c  between  adjacent  ones 

thlpeak(0)=  pi/2 

omax=  8  !  number  of  "orders"  we  consider 
600  do  601  m6  =  1,  omax 

right  =  m6*  wavelengthO/  (2*  hh) 

!  type*,  m6,  hh.  right 

if  (right  .gt.  1)  goto  602 

thlpeak(m6)=  acos(right) 

thldiff(m6)=  -thlpeak(m6)+  thlpeak(m6-l) 

620  type  621,  int(m6+.01),  thlpeak(m6) 
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235 

621 

format(’  peak’,  i4,  ’  at’,  f8.4,  ’  radians’  ) 

236 

601 

continue 

237 

602 

continue 

238 

m6max=  m6-  1 

239 

thldiff(m6max+ 1)=  thlpeak(m6max) 

240 

type  2 

241 

242 

c 

the  main  nested  loops  follow  now.  700  loop  sums  over  the 

243 

c 

reflection  peaks  (the  ’orders’);  200  loop  sums  over  pol’s 

244 

c 

(TE  or  TM  );  300  loop  over  incident  angles  around  the  peak. 

245 

700 

do  701  m7=  1,  m6max 

246 

type  2 

247 

710 

type  711 

248 

711 

format  (’  pause,  then  enter’) 

249 

accept* 

250 

720 

type  721,  int(m7  +  .01),  thlpeak(m7) 

251 

721 

format(’  peak’,  i2,  ’  at’,  f8.4,  ’  radians’  ) 

252 

c 

The  ’integrated  reflectivity"  is  found  by  summing  the 

253 

c 

reflectivity  from  thmin  to  thmax,  each  of  which  are  located  1/3 

254 

c 

of  the  way  to  the  next  peak 

255 

thcent=  thlpeak(m7) 

256 

thmin  =  thcent-  thldiff(m7  + 1)/3 

257 

thmax  =  thcent  +  thldiff(m7  )/3 

258 

thstep=  .001 

259 

44 

type  45,  thmin,  thcent,  thmax,  thstep 

260 

45 

format  (’  min,  center,  max,  step  of  theta  =’,  4f8.4) 

261 

print  2 

262 

263 

200 

do  201  pol—  1,  2 

264 

205 

if  (pol  .It.  1.5)  goto  206 

265 

type  207 

266 

207 

format  (’  TM  polarization’ ) 

267 

goto  208 

268 

206 

continue 

269 

type  209 

270 

209 

format  (’  TE  polarization’  ) 

271 

208 

continue 

272 

42 

type  43 

273 

43 

format  (’  thetal  abs{refl}  [col2j**2’) 

274 

j=-> 

275 

sum=0 

276 

300 

do  301  thetal  =  thmin,  thmax,  thstep 

277 

j=  j  +  1 

278 

theta(l)=  thetal 

279 

400 

do  401  i  =  1,  dim 

280 

sinthe(i)=  (  ninx(l)/ninx(i)  )*  sin(theta(l)) 

281 

if  (  Real(sinthe(i))  .gt.  1  )  goto  9000 

282 

!  compute  costhe  from  sinthe; 

283 

costhe(i)=  (  1-  sinthe(i)**2  )  **(.5) 

284 

p(i)=  costhe(i)*  ninx(i)/  mu(i) 

285 

!  TE  or  TM  polarization 

286 

220 

if  (pol  .It.  1.5)  goto  221 
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287 

p(i)=  costhe(i)*  mu(i)/  ninx(i) 

288 

221 

continue 

289 

beta(i)=  costhe(i)*  ninx(i)*  h(i)*  kO 

290 

401 

continue 

291 

pl=  P(l) 

292 

plast=  p(subno+2) 

293 

[plast  is  the  p  for  the  last  layer 

294 

cb2=  cos(beta(2)) 

295 

cb3=  cos(beta(3)) 

296 

sb2=  sin(beta(2)) 

297 

sb3=  sin(beta(3)) 

298 

!  insert  the  new  calculation  of  laml  and  lam2  here 

299 

jjmax=subno+  1 

300 

84 

do  85  jj  =  2,  jymax 

301 

!  define  the  transfer  matrix  for  layer  jj 

302 

tran(jj,  1,1)=  cos(beta(jj)) 

303 

tran(jj,  2,2)=  cos(beta(jj)) 

304 

tran(ij  ,1,2)=  -ii*  sin(beta(ij))/ p(ij) 

305 

tran(jj,  2,1)=  -ii*  sin(beta(jj))*  p(jj) 

306 

184 

Sprint  185,  jj 

307 

185 

[formate  tran[\  f3.0,  ’]’  ) 

308 

284 

[type  285,  tran(jj,l,l),  tran(jj,l,2) 

309 

[type  285,  tran(jj,2,l),  tran(jj,2,2) 

310 

285 

[format  (2el2.2,  ’  ’  ,2el2.2) 

311 

[type  2 

312 

85 

continue 

313 

!  type  2 

314 

!  the  product  matrix  for  the  1st  layer  (layer  2) 

315 

prod(2,  1,1)=  tran(2,  1,1) 

316 

prod(2,  1,2)=  tran(2,  1,2) 

317 

prod(2,  2,1)=  tran(2,  2,1) 

318 

prod(2,  2,2)=  tran(2,  2,2) 

319 

87 

do  88  jj  =  3,  ijmax 

320 

!  the  product  matrices  for  the  other  layers 

321 

prod(jj,l,l)=  prod(jj-l,l,l)*  tran(jy,l,l)  + 

322 

2  prod(jj-l,l,2)*tran(ij,2,l) 

323 

prod(jy,l,2)=  prod(jj-l,l,l)*  tran(jj,l,2)  + 

324 

2  prod(ij-l,l,2)*tran(jj,2,2) 

325 

prod(jj,2,l)=  prod(jj- 1,2,1)*  tran(jj,l,l)  + 

326 

2  prod(ij-l,2,2)*  tran(ij,2,l) 

327 

prod(ij,2,2)=  prod([j-l,2,l)*tran(ij,l,2)  + 

328 

2  prod(jj-l,2,2)*tran(ij,2,2) 

329 

1384 

print  385,  jj 

330 

!385 

formate  prod[\  f3.0,  ’]’  ) 

331 

!484 

type  485,  prod(jj,l,l),  prod(jj,l,2) 

332 

j 

type  285,  prod(ij,2,l),  prod(jj,2,2) 

333 

*485 

format  (2el2.2,  ’  ’  ,2el2.2) 

334 

i 

type  2 
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335 

88 

continue 

336 

!stop 

337 

a=  prod  (jj  max,  1,1) 

338 

b=  prod (jj max,  2,2) 

339 

c-  prod  (jj  max,  1,2) 

340 

d=  prodQjmax,  2,1) 

341 

342 

delta  =  b-a 

343 

discri=  (  delta**2  +  4*c*d  )**.5 

344 

f=  delta+  discri 

345 

laml  =  (a+b+discri)/ 2 

346 

lam2=  (a+b-discri)/  2 

347 

abslaml  —  abs(laml) 

348 

abslam2=  abs(lam2) 

349 

! 

type*,  laml,  lam2,  laml*lam2 

350 

type  2 

351 

c 

now  follow  two  evaluations  of  refl. 

352 

c 

490  determines  which  is  used;  the  other, 

353 

c 

though  also  exact,  might  lead  to  overflow. 

354 

psil=  2*d  +  plast*f 

355 

psi2=  f-  2*plast*c 

356 

psi3pl=  2*c*pl  +  f 

357 

psi3mi=  2*c*pl-  f 

358 

psi4pl=  f*pl  +  2*d 

359 

psi4mi=  f*pl-  2*d 

360 

361 

490 

if  (  abslaml  .It.  abslam2)  goto  491 

362 

termpl=  psi2*  psi4pl/  psil 

363 

termmi=  psi2*  psi4mi/  psil 

364 

top—  psi3mi+  termpl*lam2**(2*layno) 

365 

bot=  psi3pl+  termmi*lam2**(2*layno) 

366 

refl=  top/  bot 

367 

goto  351 

368 

491 

continue 

369 

termpl  =  psil*  psi3pl/  psi2 

370 

termmi=  psil*  psi3mi/  psi2 

371 

top=  psi4pl+  termmi*laml**(2*layno) 

372 

bot=  psi4mi+  termpl*laml**(2*layno) 

373 

refl=  top/  bot 

374 

375 

351 

continue 

376 

c 

no  longer  needed  avrefl(j)=  abs(refl) 

377 

rr=  abs(refl)**2 

378 

!  print  only  values  above’floor" 

379 

floor  =  .0 

380 

70 

if  (abs(refl)  .It.  floor)  goto  71 

381 

c 

if  (mod(j,10.)  .gt .  0)  goto  71 

382 

c 

preceding  line  cuts  down  on  the  printout 

383 

c66 

print  67,  thetal,  abs(refl),  rr 

384 

67 

format  (2f8.4,  f9.5) 

385 

71 

continue 

386 

sum=  sum+  rr 

19 


387 

301 

continue 

388 

mnteg=  thstep*  sum 

389 

50 

type  51,  rrinteg 

390 

51 

formate  integrated  reflectivity  =  el2.3  ) 

391 

302 

goto  303 

392 

9000 

type  9001,  thetal 

393 

9001 

format  (’  no  penetration  for  thetal  greater  than’,f8.4) 

394 

303 

type  2 

395 

201 

continue 

396 

701 

continue 

397 

!90 

type  91, floor 

398 

91 

format  (’  floor  =  ’,  f3.2) 

399 

400 

9999 

end 

401 

402 

403 
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Figure  2.  Beam  entering  a  double  layer. 
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Figure  3.  Beam  entering  N  double  layers 
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